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Introduction 



This is the second installment of a report-pair concern- 
ing implementation of tensor product factoring of coeffi- 
cient matrices in applications of the finite element method 
to numerical weather prediction. It was noted in Part 1 
(Ref. 1) that these techniques were introduced in numerical 
weather prediction by Staniforth and Mitchell (Ref. 2). 
Discussed in Part 1 are applications in which the "mass" 
matrix for a grid such as that shown in Fig. 1 is factored 
as the tensor product of two matrices. 
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1. Node numbering and spacing. 



One of these matrices (MA) depends solely upon the nodal 
spacing in the east-west direction (a^) and the other (MB) 
depends only on the north-south spacing (b^)* We began with 
the set of simultaneous linear equations 

M w = v , < 1> 

where M (the "mass" matrix) is symmetric, ne x ne , and w and 
v are column vectors of height ne . M and v are input quan- 
tities and w is sought. The tensor product representation 
of M is 



M = MB * MA, <2> 

where MB and MA are tridiagonal, symmetric matrices, e x e 
and n x n, respectively. (The tensor product and matrices 
MA and MB are defined in Appendix A. ) This representation 
allowed <1> to be rewritten as 



MA W MB = V, <3> 

where W is n x e and the successive columns are subvectors 
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V is also n x e 



of w corresponding to the rows of Fig. 1. 
and similarly derived from v. Boundary conditions consid- 
ered were a cyclic condition in the east-west direction and 
either homogeneous Neumann conditions (normal derivative 
zero) or nonhomogeneous Dirichlet conditions (specified 
nonzero values) on the northern and southern edges. 

The present report discards the cyclic east -west boundary 
condition and deals with two cases: 

(1) Solutions of <3> with nonhomogeneous Dirichlet con- 
ditions on all four edges ; 

(2) Solution of Poisson's equation for the same region 
with nonhomogeneous Dirichlet conditions on all four 
edges . 



Mass Matrix Dirichlet Boundary Conditions 



Effects of the Dirichlet boundary conditions on the solu- 
tion process are most readily understood by considering the 
following partitioned form of <1>: 



In <4> the w 
boundary values 



[Hh Hi® = [$ 

vector has been rearranged so that 



are in the subvector 



Wi 



and 



<4> 

all of the 
the interior 



("center”) values are in w c . A similar reordering has been 
applied to v and M. If the boundary values of w are pre- 
scribed, then w^ is known and only w c remains to be found. 
Expanding the lower partition of <4> and placing the known 
terms on the right gives 

M 22 W C = v c ‘ i w b » <3> 



or, letting v ’ = v c - M 2 iWj 3 , we have 

M 22 w c = v c ’ . 



<5 ' > 



We consider now how the strategy just described can be 
applied when the tensor product resolution of M has been 
used to convert <1> into <3>. In the matrix W the pre- 
scribed boundary values occupy the first and last columns 
and the top and bottom rows. Denote this border matrix, 
including an (n-2) x (e-2) null matrix inside, by WB . 
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Calculate 



VB = MA WB MB 



<b> 



and now form 



V' = V - VB. 



<7> 



Now define a set of submatrices MAI, MB1, Wl, and VI 
obtained from MA, MB, W, and V', respectively, by removing 
the first and last columns and the top and bottom rows. The 
reduced problem becomes 

MAI Wl MB1 = VI <8> 

As described in Ref. 1, <8> may be solved by standard Gaus- 
sian elimination procedures. A computer program (GAUSS4) 
which carries out these calculations is listed in Appendix 
B. The subroutines of GAUSS4 are designed for substitution 
in the program devised by Hinsman (Ref. 5). 



Poisson's Equation Dirichlet Boundary Conditions 

As noted above, Staniforth .and Mitchell (Ref. 2) appear 
to have been first in applying the tensor product resolution 
to Poisson’s equation in a numerical weather prediction 
problem using the finite element method. Additional detail 
is given in earlier papers by Dorr (Ref. 3) and by Lynch, 
Rice, and Thomas (Ref. 4). 

Finite element discretization of Poisson's equation for 
the region of Fig. 1 results in a set of simultaneous linear 
equations which may be written in matrix form as 

K w = v, <9> 

where vectors v and w are, respectively, given and unknown. 
As for <1>, each has length ne and the coefficient matrix K 
is ne x ne , symmetric, sparse, and block tridiagonal. K is 
called the "stiffness" matrix in finite element parlance. 

It is easily shown that K is expressible as the sum of 
two tensor products as follows: 

K = MB * SA + SB * MA. <10> 

The new matrices SA and SB are symmetric, tridiagonal and 
depend only on the a., and b., respectively. Explicit formu- 
las for SA and SB are given in Appendix A. 
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Using the definition of the tensor product and again con- 
verting the vectors w and v into the n x e rectangular 
matrices W and V, <9> may be written as 

SA W MB + MA W SB = V. <11> 

Before solving <11> we must first take account of the Diri- 
chlet boundary conditions on the four edges of the region. 
As in solving <3>, the given boundary values are in the 
first and last columns and top and bottom rows of W. As 
before, we let WB be an n x e matrix containing the given 
boundary values , together with zeros at locations corre- 
sponding to interior nodes. Calculate 

VB = SA WB MB + MA WB SB, <12> 

and then form 



V’ = V - VB. <7> 

The remaining step again parallels that used when applying 
the Dirichlet boundary conditions to <3>. Specifically, we 
introduce submatrices MAI, MB1, SAl , SB1, W1 , and VI 
obtained from MA, MB, SA, SB, W, and V', respectively, by 

0 i 

removing the first and last columns and the top and bottom 
rows. The reduced problem becomes 

SAl W1 MBl + MAI W1 SBl = VI. <13> 



To solve <13> we first need the complete solution of the 
eigenproblem 

SBl p.j = A. j MBl p.j, <14> 

where p^ is the ith eigenvector and A • is the corresponding 
eigenvalue. We write the complete solution in the form 

SBl P = MBl P A, < 14 ' > 

where P is the (e-2) x (e-2) modal matrix whose columns are 
the p.. and A is the (diagonal) spectral matrix whose ele- 
ments are the Aj . We specify that the modal matrix is nor- 
malized so that 

PT MBl P = I, < 15> 

where I is the identity matrix of order e-2 and PT is the 
transpose of P. If both sides of <13> are postmultiplied by 
P and < 14 ' > is used to replace SBl P, <13> becomes 

SAl W1 MBl P + MAI W1 MBl PA = VI P. <16> 

Let X = W1 MBl P and U = VI P , then <16> is equivalent to 
(SAl + A .j MAI) x.. = u^. , i = 1, e-2, <17> 
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where x.. and u. are, respectively, the ith columns of X and 
U. Since the coefficient matrix in <17> is tridiagonal, the 
Gaussian elimination process, i.e., factoring, forward 
reduction, and back- subs titution , is computationally econom- 
ical. The final step consists of a matrix multiplication to 
obtain 

W1 = X PT. <18> 

Since Wl contains the w values at all interior nodes and the 
boundary values were known in advance, the solution is com- 
plete. A FORTRAN program (GAUSS5) which implements the ten- 
sor product solution for Poisson's equation is given in 
Appendix C. 



Operation Counts and Storage Requirements Poisson ' s Equa - 
tion 

In Ref. 1 comparisons of floating point operation counts 
and storage requirements were made for solutions of <1>. 
Substitution of the boundary conditions considered here in 
place of those considered in Ref. 1 has a negligible effect 

t , 

on both operation counts and storage requirements. Accord- 
ingly, no further comparison is given here for solutions of 
< 1 > . 

Solution of Poisson's equation (<9>) using the tensor 
product resolution <10> of K is more costly in terms of com- 
putation and storage than the previously studied applica- 
tions to <1>. In Table 1 the number of floating point oper- 
ations and the required number of coefficient matrix storage 
locations are compared for three different solution methods. 
These are SOR (successive over- relaxation) , SKY (skyline 
storage and Gauss elimination), and TENSOR (the scheme 
described above). A floating point operation is defined to 
be one multiplication (or division) plus one addition (or 
subtraction). The exact operation counts would be polynomi- 
als in n and e. Only the highest degree terms are given in 
the table. Since it is not possible to predict the number 
of iterations per solution using SOR, the operation count 
given for that algorithm is for a single iteration . In 
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Table 1 a storage location corresponds to 8 bytes. For the 
comparison it is assumed that each floating point number 
requires 8 bytes of storage and an integer requires 4 bytes. 
The storage requirement given for SOR is based on the com- 
pact storage scheme described by Franke and Salinas 
(Ref. 6). 



TABLE 1. Operation Counts and Storage Requirements. 



ALGORITHM 


NUMBER OF OPERATIONS 


NUMBER OF STORAGE LOCATIONS 




PER SOLUTION 


FOR COEFFICIENT MATRICES 


SOR 


10 en (1) 


13 en 


SKY 


2 en 2 


en 2 


TENSOR 


2 en 2 


e 2 



Note: 1. Number of operations per iteration . 



It is perhaps surprising to note that the number of oper- 
ations for TENSOR is no fewer than for SKY. Turning atten- 
tion to storage requirements reveals that for a large prob- 
lem (e = n = 100, say) the SKY storage requirement for the 
stiffness matrix is 8 megabytes, compared with 1 megabyte 
for SOR and 80 kilobytes for TENSOR. It is this comparison 
which is the compelling reason for preferring TENSOR. It is 
acknowledged that there is overhead associated with the one- 
time solution of the eigenvalue problem <14>, but the tri- 
diagonal form of matrices SBl and MB1 makes the amount of 
computation comparable with that required for a single solu- 
tion of Poisson's equation. Since two solutions of Pois- 
son's equation are required at each time step, the overhead 
is clearly negligible. 

It is not feasible to make a definitive comparison 
between the number of operations required for SOR and those 
required for the other two algorithms. If the number of 
iterations is less than 0.2 e, then SOR will be more econom- 
ical and the storage tradeoff would need to be weighed. 
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Conclusions 



It has been demonstrated that Dirichlet boundary condi- 
tions on all edges of the region are easily incorporated in 

solution processes which use tensor product resolution of 
the coefficient matrix. For very large problems the tensor 
product algorithm uses much less core storage than alterna- 
tive choices . The computational expense of a solution to 
Poisson's equation is substantially the same for Gaussian 
elimination and for the tensor product scheme. It is 
expected that successive over-relaxation is almost always 
more expensive. 
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APPENDIX A - TENSOR PRODUCT AND MATRIX DEFINITIONS 



Tensor Product The tensor product of matrices C and D 

may be represented in block partition form as 





■cn D 


C12 D 


c 1 3 D“ 


C *0 = 


C21 D 


C22 D 


C23 D 




_c 3 i D 


c 3 2 D 


c 3 3 P_ 



where the c-jj are the elements of C. Note that, if C and D 
have dimensions r x s and t x u, respectively, the tensor 
product has dimensions rt x su . 



Definitions for matrices MA and SA are given below. The 
corresponding expressions for MB and SB may be obtained by 
substituting "b" for "a" throughout and replacing n by e . 
(Symbols a-j and b -j are defined in Fig. 1.) 



(n=4) 



SA = 
(n= 4 ) 



2a x 


ai 


0 


0 


^1 


2(a!+a 2 ) 


a 2 


0 


0 


a 2 


2(a 2 +a 3 ) 


a 3 


0 


0 


a 3 


2a 3 



’ I 
ai 


1 

‘ a x 


0 


0 




I + I 




0 


ai 


a v a 


2 a 2 


0 


_1 


1+1 . 


1 


a 2 


a 2 a 3 


a 3 


0 


0 


a 3 


l 

a 3 
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APPENDIX B 



PROGRAM to SOLVE M w = v with DIRICHLET BOUNDARY CONDITIONS 
Listing: GAUSS4 FORTRAN 



C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 



MAIN PROGRAM MASS MATRIX USING TENSOR PRODUCT RESOLUTION 



THIS PROGRAM IS DESIGNED TO TEST THE SCHEME (TENSOR) 

A TENSOR PRODUCT 

ORDER TO SOLVE THE SYSTEM OF EQUATIONS M w = v . IN 



WHICH RESOLVES THE MASS MATRIX INTO A 



IN 



1000 



THIS PROGRAM THERE ARE DIRICHLET BOUNDARY CONDITIONS ON 
ALL 4 EDGES OF THE REGION. THE PRESCRIBED BOUNDARY 
VALUES ARE GIVEN IN THE CORRESPONDING LOCATIONS IN V. 

THE SUBROUTINES MAY BE INSERTED IN THE PROGRAM DEVISED 
BY HINSMAN. 

IMPLICIT REAL-8 (A-H.O-Z) 

COMMON / CM1A/ NLAT , NLONG 
COMMON / CM8 / A ( Z 1 ) .B(Z1) 

COMMON AG(ZB ) ,BG(ZC) , GAD ( ZK ) , GBD ( ZL ) ,MA(ZM) ,MB(ZN) 
DIMENSION V(ZP) 

READ(5 , *) NLONG, NLAT 
LATX=NLAT+ 1 
WRITE (6, 1000) 

FORMAT ( / / , ' MASS MATRIX - TENSOR PRODUCT RESOLUTION' 



1TE ( 6 , 1001) NLONG, NLAT 

READ (5 . - )A.B 
WRITE(6 ,506)A 
WRITE(6,503)B 

FORMAT (/ , ' B: ’,(24F3.0)) 

FORMAT (/,' A: ' ,(24F3.0)) 



503 
500 

1001 FORMAT NLONG = \I3,' " NLAT =’,13,/) 

C CONSTRUCT FACTORS, GAD AND GBD, OF MASS MATRIX 
CALL AMTRX3 
WRITE (6. 501) AG 
FORMAT?// AG: 

WRITE (6 ,504 )BG 
FORMAT?/,’ BG: 

WRITE ? 6 , 1002 )GAD 
FORMAT? / , ' GAD ' , / , ( 3X , 6F7 . 3 ) ) 
WRITE(6.i004)GBD 
FORMAT?;/. :_GBD' , / , ( 3X , 6F7 . 3 ) ) 

WRITE ( 6 , 10 



501 

504 

1002 

1004 



\(12F4.1)) 
, ( 12F4 . 1 ) ) 



C 

4 



, , _ 03 )MA 

WRITE (6 , 1006 JMB 
K= (NLAT- 1) -NLONG 
L= NLAT -NLONG 
READ (5 , - ) V 
WRITE? 6 , 5 10 ) V 

CORRECT RIGHT-HAND SIDE FOR DIRICHLET CONDITION 
LONGM=NLONG- 1 
DO 2 J=2 , LONGM 

V ( J + NLONG ) = V ( J + NLONG ) - ( GAD? 2* J - 1 ) *V ( J - 1 ) + GAD ( 2* J - 2 ) 

1 - V ( J ) + GAD ( 2 - J + 1 ) - V ( J+ 1 ) ) -GBD ( 3 ) 

V(K + J)=V(K+J3- (GADC2 * J- 1 ) -V (L+ J- 1 ) +GAD (2* J- 2 ) *V(L+ J ) 

1 + GAD? 2 - J + 1 ) *V (L+ J + 1 ) ) -'GBD ( 2-LATX - 1 ) 

CU=GBD(3) 

CL= GBD (2-LATX- 1) 

GBD?3 ) =0 . 

GBD f 2-LATX- 1)=0 . 

DO 3 J = 2 , NLAT 

V ? ( J - 1 ) -I?LONG + 2 ) = V ( ( J - 1 ) -NLONG + 2 ) - ( GBD ( 2* J - 1 V *V ( ( J - 2 ) 
1- NLONG + 1) + GBD ( 2 - J - 2 J -V ( ( J - 1 ) -NLONG + 1 ) + GBD ( 2* J + 1 ) 

2 - V? J -NLONG + 1 ) j '-'GAD ( 3 ) 

V ( J - NLONG - 1 ) =V (J -NLONG - 1 ) - ( GBD ( 2 - J - 1 ) - V ( ( J - 1 ) -NLONG ) 

1 + GBD ( 2* J - 2 ) -V ( J*NLONG ) + GBD ( 2* J+ 1 ) *V ( ( J+ 1 ) -NLONG )) 

2 * GAD ( 2* NLONG - I ) 

GBD(3)=CU 

GBD ( 2-LATX- 1)=CL 

WRITE (6, 5 10 )V 

PERFORM LDLT FACTORING OF GAD AND GBD 
CALL FACT 1( GAD, NLONG) 

CALL FACT 1 ( GBD , LATX ) 

WRITE (6, 1002 )GAD 
WRITE (6, 1004 )GBD 
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c 

c 



c 

c 



6 

510 

1003 

1006 



PERFORM FORWARD REDUCTION AND BACK- SUBSTITUTION USING 
FACTORS OF GAD 

CALL BACKA 1 ( GAD , V ) 

WRITE (6, 5 10 JV 

PERFORM FORWARD REDUCTION AND BACK- SUBSTITUTION USING 
FACTORS OF GBD 

CALL BACKB 1 ( GBD , V ) 

WRITE?6 .510 )V 



FORMAT () ,' V: ' , 5F8 . 2 . / , f 4X , 5F8 . 2 ) ) 

format /;; MA:;,2X,$6i3' 



C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 



C 

C 

C 



FORMAT (/ , 
STOP 
END 



mb: ' ; 2x; 3613) 



SUBROUTINE FACT1(A,NN) 



SUBROUTINE FACTl PERFORMS L*D*LT FACTORING ON A SUBMATRIX 
OF A SYMMETRIC TRIDIAGONAL MATRIX STORED IN SKYLINE FORM. 
THE SUBMATRIX IS FORMED BY OMITTING THE FIRST AND LAST 
COLUMNS AND ROWS OF THE INPUT MATRIX. 

, - - INPUT VARIABLES 

A(NWK) = INPUT MATRIX STORED IN COMPACTED FORM 
NN = NUMBER OF COLUMNS (OR ROWS) IN INPUT MATRIX 

NWK = NUMBER OF ELEMENTS BELOW SKYLINE (2*NN - 1) 

- - OUTPUT - - 

A (NWK) = D AND L - FACTORS OF INPUT SUBMATRIX 



IMPLICIT REAL*8 (A-H, 0- Z ) 

DIMENSION A( 1) 

PERFORM L*D*LT FACTORIZATION OF STIFFNESS MATRIX 

LONGMM=NN- 2 
A( 3 ) = 0 . 

DO 50 J=2 ,LONGMM 



120 

50 

2000 



TEMP = A ( 2 * J + 1 ) / A ( 2 * ( J - 1 ) ) 
A ( 2*J)=A( 2*J j - T&MP -A ( 2 - J + 
IfTaT 2' v JJ 1120, 120,50 
WRITE (IOUT, 2000 )N,A(KN) 



1 ) 



STOP 

A ? 2* J + 1 ) = TEMP 

FORMAT (/ / , ' STOP - MATRIX NOT POSITIVE DEFINITE’ //, 
1' NONPOSITIVE PIVOT FOR EQUATION ',14,//,' PIVOT = ' 
2E20.12) 

RETURN 
END 






U ^ ^ ^ ^ ^ 



C 

C 

c 



c 

c 

c 



20 

C 

C 

c 

40 

C 

C 

C 



4 % 4 \ 4 % 4 \ 4 \ S* 4 % 4 \ 4 % 4 % 4 \ 4 \ 4 \ 4 % 4 % 4 \ 4 % 4 % 4 % 4 % 4 % 4 \ 4 % 4 \ 4 \ 4 \ 4 % 4 \ 4 % 4 \ 4 \ 4 \ 4 \ 4 \ 4 % 4 % 4 \ 4 % 4 \ 4 % 4 % 4 % 4 \ 4 \ 4 \ 4 % 4 \ 4 \ 4 \ 4 \ 4 \ 4 % 4 % < 

THIS SUBROUTINE PERFORMS THE FORWARD REDUCTION AND BACK- 
SUBSTITUTION USING THE FACTORS OF GAD 

IMPLICIT REAL*8 (A-H.O-Z) 

COMMON 7 CM 1A/ NLA T . NLONG 
DIMENSION A( 1 ) ,V(1) 

DEFINE LIMITS FOR DO- LOOPS 

NTM=NLAT- 1 
LONGM= NLONG -1 
LONGMM= NLONG - 2 

REDUCE RIGHT-HAND- SIDE LOAD VECTOR 

DO 100 K= 1 , NTM 
DO 20 J=3,LONGM 

V (K*NLONG + J ) = V (K* NLONG + J ) - V (K -NLONG + J-1)*A(2*J-1) 
DIVIDE BY DIAGONAL ELEMENTS 
DO 40 J= 1 , LONGMM 

V (K* NLONG + J + 1 ) = V ( K -NLONG + J + 1 ) / A ( 2 * J ) 

BACK- SUBSTITUTE 
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DO 60 J=3 , LONGM 
L= (K + lJ-'NLONG- J + 1 
M= 2* ( NLONG - J ) + 3 

60 v(l)=vTl)-v(L+i)*a(m) 

100 CONTINUE 

RETURN 

END 



C 

c 

c 

c 



c 

c 

c 



c 

c 

c 



20 

c 

c 

c 

40 

C 

C 

C 

60 

100 



C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 



SUBROUTINE BACKB1 (A , V ) 

^ ^ ^ # ^ «k*# %*<• %*# %*<• %*# K*# 

/« /% /% 4 % /« /« #« #% /« /% /% /% 4 % /% 4 % #% #% 4 % /% /% #% 4 % /« /« /% #« /% #% /« /% /% /\ /\ /% /% #« 

THIS SUBROUTINE PERFORMS THE FORWARD REDUCTION AND BACK- 
SUBSTITUTION USING THE FACTORS OF GBD . 

IMPLICIT REAL*8 ( A-H. 0- Z ) 

COMMON/ CMlA/NLAt .NLONG 
DIMENSION A(1),V(1) 

DEFINE NEEDED INDEX VARIABLES 

LATX=NLAT+ 1 
LONGM=NLONG- 1 

REDUCE RIGHT-HAND- SIDE LOAD VECTOR 

DO 100 K=2, LONGM 
DO 20 J = 3 , NLAT 

V (K+ ( J- 1)“ NLONG )=V(K+(J-1) *NLONG ) - V (K+ ( J - 2 ) *NLONG ) 

1 - A ( 2 - J - 1 ) 

DIVIDE BY DIAGONAL ELEMENTS 
DO 40 J = 2 , NLAT 

V ( K + ( J - 1 ) -NLONG ) = V ( K + ( J - 1 ) -NLONG ) / A ( 2 * J - 2 ) 

BACK- SUBSTITUTE 

DO 60 J = 3 , NLAT 

V (K + ( LATX - J T -'NLONG )=V (K+ ( LATX - J ) -'NLONG ) 

1 - V ( K+ T 1 + LATX - J ) -NLONG ) - A ( 2 - ( LATX - J) + 3 ) 

CONTINUE 
RETURN 

* <t ** / * / * /» /« #* /\ /* <« » » /» < * / » # » »» *\ /« ** '*/*/*'*«*<*/*** **'*'*<* **/w* /*'*/* * ******* /*<****» '* /» /* /\ /* /« < * * * / * «» /* 

% 4 % 4 % 4 % 4 % #% 4 % 4 \ 4 \ 4 \ 4 % 4 \ 4 \ 4 > 4 % 4 \ 4 % 4 % 4 \ 4 % 4 % 4 \ 4 % 4 % 4 % 4 % * 4 % * 4 % * 4 % 4 % 4 % * 4 % * 4 % * 4 % 4 % 4 % * 4 % * 4 \ 4 % 4 % 4 % 4 \ 4 % 4 % 4 \ A 4 \ 4 \ 4 % 4 % 4 % 4 % 4 % 

THIS SUBROUTINE FORMS THE MASS MATRIX IN THE FORM OF A 
TENSOR PRODUCT OF THE GBD MATRIX AND THE GAD MATRIX. 

THE FIRST OF THESE IS NLAT + 1 BY NLAT + 1, SYMMETRIC, 

AND TRIDIAGONAL. THE SECOND IS NLONG BY NLONG, SYMMET- 
RIC, AND TRIDIAGONAL. NOTE THAT THERE ISNO CYCLIC 
BOUNDARY CONDITION IN THE EAST-WEST DIRECTION. BOTH GBD 
AND GAD ARE STORED IN SKYLINE VECTOR FORM ?UPPER TRIANGLE 
WITH SPACE FOR FILL-IN). INTEGER ADDRESS VECTORS MB AND 
MA ARE ALSO GENERATED. 



C 

C 

c 



2 

C 



IMPLICIT REAL*8 (A-H ,0- Z ) 
COMMON/ CM 1A/ NLAT , NLONG 
COMMON/ CM8 /A( Z1 ) .bTzI) 
COMMON AG(ZbLbG(Z 
DIMENSION BG ( NLAT ‘ 

1GAD ( 3 -NLONG- J ) ,MA 

LATX = NLAT + 1 
LONGM= NLONG -1 

FIND BG = (ELEMENT HEIGHT)/ 6 
LONGM= NLONG -1 
DO 2 J= 1 , NLAT 
BG ( J ) =B ( i+LONGM* ( J- 1) ) / 3 . 
GENERATE GBD 

GBD(1)=2.*BG(1) 

DO 4 J=2 , NLAT 
K=2- (J- 1 ) 

GBD (K) =2 . * (BG ( J- 1 ) +BG ( J ) ) 
GBD(K+1)=BG(J-1) 

GBD ( 2 -NLAT ) =2 . -BG ( NLAT ) 

13 



B(Z1 ) 

ZC ) , GAD (ZK) , GBD ( ZL ) , MA ( ZM ) ,MB(ZN) 
) , AG (NLONG) , GBD (2*NiAT-l) , 
(NLONG+1) , M^ ( NLAT + 2 ) 



4 



c 

10 

c 

12 

C 

16 

18 



GBD ( 2*NLAT + 1 ) = BG ( NLAT j 
FIND AG = (ELEMENT WIDTH) /6. 
DO 10 J= 1 , LONGM 
AG ( J ) =A( J ) / 3 . 

GENERATE GAD 

GAD ( 1 ) = 2 . *AG ( 1 ) 

~~ 12 J = 2, LONGM 



DO 12 J=2, LONGM 

K=2*?J-1) 

GAD(K)=2.*(AG(J- 



GAD 

GAD 

GAD 



1 ) + AG ( J ) ) 



K+1]=AG(J-1) 

2 “LONGM ) = 2 *AG (LONGM ' 
' 2*LONGM+ 1 ) = AG (LONGM ' 
GENERATE DIRECTORY VECTOR^ 
MB ( 1 ) = 1 

DO 16 J = 1 , NLAT 
MB(J+lI=2“J 

NLAT + 2 ) = 2* ( NLAT + 1 ) 



MB I 
MAI 
DO 



= 1 

J = 2 ,NLONG 



MA( JJ =2* ( J- 1 ) 
MA1NLONG+ 1 ) =2*NLONG 
RETURN 
END 
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APPENDIX C 

PROGRAM - POISSON'S EOUATION with DIRICHLET BOUNDARY 

CONDITIONS 



Listing : GAUSS5 FORTRAN 

MAIN PROGRAM 



STIFFNESS 

RESOLUTION 



MATRIX USING TENSOR 



THIS PROGRAM IS DESIGNED TO TEST THE SCHEME WHICH 
RESOLVES THE STIFFNESS MATRIX INTO A SUM OF TWO TENSOR 
PRODUCTS IN ORDER TO SOLVE THE SYSTEM OF EQUATIONS 
K W = V. THERE ARE DIRICHLET BOUNDARYCONDITIONS ON 
ALL 4 EDGES OF THE REGION. THE PRESCRIBED BOUNDARY 
VALUES ARE GIVEN IN THE CORRESPONDING LOCATIONS IN V. 
THE SUBROUTINES MAY BE INSERTED IN THE PROGRAM DEVISED 
BY HINSMAN. 

IMPLICIT REAL*8 (A-H.O-Z) 

COMMON/ CM1A/NLAT,NL0NG 



Z1KB(Z1) 

,BG(ZC) ,GA1(ZK) , SAl (ZK) .GBl(ZL) ,SB1(ZL) 

p) ,wi(zd) , p(zr) d(zs) ,u?zt) 



COMMON/ CM8/ A 
COMMON AG (ZB 

DIMENSION V(__ . 

READ ( 5 , *)NLON6, NLAT 
LATX=NLAT+ 1 
WRITE (6 ,1000) 

1000 FORMATT//,’ STIFFNESS MATRIX - TENSOR PRODUCT 
1RESOLUTION' /) 

WRITE? 6 , 100 1 )NLONG , NLAT 
READ (5 .*)A,B 
WRITE (6 ,500) A 
WRITEf 6 ,503)B 
FORMAT?/,' B: ' , (24F3 

FORMAT?/,' A: ' , ( 24F3 

NLONG = ,13, 

FACTORS, GAl, GB1, 



503 
500 
1001 
C 
C 



501 

504 

1002 

1012 

1004 



FORMAT ( 
CONSTRUCT 
MATRIX 

CALL AMTRX4 
WRITE (6 , 50 1 ) AG 

TTniDVr A r T f 7 T A 



ojj 



SAl 



NLAT =’ 13./) 

AND SA2 OF STIFFNESS 



FORMAT t / , a 
WRITE?6 ,504)BG 
FORMAT?/,’ BG: 

WRITE? 6 , 1002 ) GAl 
FORMAT?/,' GAl',/ 

WRITE?6 , i012)SAl 
FORMAT?/ ' SAl',/ 

WRITE? 6 , 1004 )GB1 
FORMAT?//.' GB1 ' , / , ( 3X , 6F7 
WRITE? 6 , l6 14 ) SB1 

-r FORMAT?//,' SB1’ / , (3X, 6F7 . 3) ) 
LOAD BORDER VECTOR wi 
READ ( 5 , ■ * V V 
CALL BORDER (W1,V) 

WRITE(6 , 510;V 
L1=4"NL6nG 
L2=L1+ 1 
L3=L1+4*LATX 



AG: 



1014 
C 



, ( 12F4 . 1 ) ) 

, (12F4.1)) 
(3X, 6F7 . 3 ) ) 

( 3X , 6F7 . 3 ) ) 

3)) 



520 

521 
510 

C 

C 

C 

C 

C 



WRITEf 6 ,520) (W1(L) ,L=1,L1) 
WRITE?6 ,521) (Wl(L) ,L=Li.L3) 
FORMAT (// , ’ W1 V ,/ (3X,^F8.2)) 
FORMAT (3X5F8. 2) 



FORMAT ( 
CALL MU 
WRITE 
WRITE 
WRITE 
CALL 
WRITE 
WRITE 
CALL 
WRITE 
WRITE 
WRITE 
READf 



Li i 



V: ’ ,/ , (4X.6F10.5)) 

a,v,dii;sBi) 



k (^ i?l) , i=i;ti) 
(W1?L) ,L=L2 ,L3 ) 

XhU 



)V 



,L= 1 , LI) 
L=Li , L3 ) 



530 



6,520 
6 521 
6*510 
ORDER 
6,520, , , 

6,521j(Wl(L, ,u-Lu.. 
uiTlJ^i,V,SAi,GBl) 

6*520 
6 j 521 

WRITtf^ Jfo)P 

FORMAT?/,' P 



n^,L lh Li) 



W1(L 



L3) 



,/ , (6X, 3F12 . 4) ) 
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READ ( 5 , * ) D 
WRITER ,531)D 
FORMAi ( / ,’ D : 



531 
C 

C FORM U = VINT*P 
C 



3F12.4) 



28 

29 

30 

532 



C 

C 

C 



38 

39 

40 



C 

C' 

C‘ 

C 

C 

C 

c 



L0NGM=NL0NG- 1 
LONGMM=NLONG- 2 
NLATM=NLAT- 1 
DO 30 L= 1 , NLATM 
JP1=(L-1) "NLATM 
KU1= f L- 1 1* (NLONG-2 ) - 1 
DO 29 K=2 , LONGM 
TEMP=0 . 

DO 28 J=l. NLATM 
JV=J*NLONG+K 
JP=JP1+J 

TEMP = TEMP + V ( JV ) *P ( JP ) 

U(KU1+K)=TEMP 
CONTINUE 
WRITEX6 ,532)U 

FORMAT (7 ’ U: ’ /, (6X.4F12.4)) 

CALL FFFf)B(U,D,GAi,^Al) 

WRITE ( 6 , 532)6 

PUT FINAL RESULT IN V 

DO 40 L=l, NLATM 
DO 39 K= 1 , LONGMM 
TEMP=0 . 

DO 38 J=l, NLATM 

TEMP = TEMP + P ( ( J - 1) "NLATM +L ) *U ( ( J - 1 ) *LONGMM+K ) 

v7l*nlong+k+i)=temp 

CONTINUE 
WRITE (6 , 5 10 ) V 
STOP 
END 



SUBROUT INE BORDER (W1 , V ) 



THIS SUBROUTINE CLEARS THE BORDER VECTOR Wl AND 
SUBSTITUTES THE BOUNDARY VALUES FROM V. 

IMPLICIT REAL'" 8 (A - H ,0 - Z ) 

COMMON/ CM lA/NLAT, NLONG 
DIMENSION Wl(ZQ) ,V(ZP) 

LATX=NLAT+ 1 

NC=4*NLONG 

NR=4'"LATX 

NB=NC+NR 

DO 4 J= 1 , NB 

W1(J)=0. 

DO 6 J= 1 .NLONG 

wi(j)=vm 

DO 8 J=l, NLONG 
L=3"NLON6+J 
K=NLAT" NLONG +J 
Wl (L)=V (K) 

DO 10 J= 1 ,LATX 
L=NC+ J 



10 



12 



K= (J-l)"NLONG+l 
‘ (L)=V (K) 

12 J=1 ,LATX 



Wl 
DO 

L=NC + 3 "LATX + J 
K=J*NLONG 

wia)=v(K) 

RETURN 
END 






c 

c 

c 

c 



SUBROUTINE AMTRX4 

+ ^ 

k /% /% #% #% #% #% «« #\ #% /% S* /% #% /% /% /% /% /% #% #% 

THIS SUBROUTINE FORMS THE MATRICES GA1 , GB1, SA1 , AND SBl 
THAT ARE FACTORS IN THE TENSOR PRODUCTS USED TO FORM THE 
COEFFICIENT MATRIX ("STIFFNESS” MATRIX) FOR THE POISSON 
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c 

c 

c 



c 

c 

c 



2 

C 



EQUATION. ALL OF THESE MATRICES ARE 
TRIDIAGONAL . 



SYMMETRIC AND 



6 

C 

10 

C 
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IMPLICIT REAL-'S ( A-H ,0- Z ) 
COMMON/ CM 1A/ NLAT , NLONG 
COMMON / CMS / A ( Z 1 > . B ( Z 1 ) 
COMMON' AG ( ZB } , BG ( ZC ) , GAl ( ZK 
DIMENSION BG( NLAT), AG (NLONG 
1GA1 ( 3 “NLONG- 3 ) 

LATX=NLAT+ 1 
LONGM= NLONG -1 

FIND BG = (ELEMENT HEIGHT) /6. 
NM= NLONG -1 
DO 2 J = 1 , NLAT 
BG(J)=B(1+NM*(J-1) )/3. 
GENERATE GBI AND 6*SB1 
GBI ( 1 ) =2 . *BG ( 1 ) 

SB 1 ( 1 J = 1 . /BG ( 1 ) 



) , SAl ( ZK ) , GBl(ZL) 
) ;GB1(2*NLAT-I) , 



(ZL) 



DO 

K=2 

GBI 

GBI 

SB1 

SB1 

GBI 

GBI 

SB1 

SB1 



4 J=2 , NLAT 
“ (J-l) 

(K)=2. *(BG(J-1)+BG(J) ) 
,K+ 1) =BG?J- 1) 

K ) = 1 . /BG ( J- 1)+1./BG(J) 

' 2 “NLAT ) = 2 . “BG (NLAT ) 

' 2 “NLAT + 1 ) = B G ( NLAT ) 
'2*NLAT)=1. /BG(NLAT) 



^ ~ ^ / 2 *NLAT + 1) =- 1 . / BG (NLAT ) 
J2=i*NLAT+l 
DO 6 J=1,J2 
SBl(J) = SBl(J)/6 . 

FIND AG = (ELEMENT WIDTH) /6. 
DO 10 J= 1 , LONGM 
AG ( J ) =A( J J / 3 . 

GENERATE GAl AND 6*SA1 
GAl ( 1 ) = 2 . “AG ( 1 ) 

SAl( 1 ) = 1 . / AG (1 ) 

DO 12 J=2, LONGM 
K=2*(J-1) 

_ (K ) =2 . * (AG ( J- 1 ) +AG ( J ) ) 
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C, , 
C“ ' 

c 

c 

c 

c 

c 

c 

C“‘ 

c 



GAl 
GAl 
SAl 
SAl 
GAl 
GAl 
SAl 
SAl 
J2 = 
DO 



K+ 1 ) =AG ( J- 1 ) x 
K ) = 1 . / AG ( J - 1 ) + 1 . / AG ( J ) 
K+1)=-1./AG(J-1) 

2 “LONGM ) = 2 “'AG ( LONGM ) 

2 “LONGM+ 1 ) = AG (LONGM) 

2 “'LONGM) = 1 . / AG(LONGM) 
2*L0NGM+1)=- SAl (2* LONGM) 
_*NL0NG-1 
14 J=1,J2 



SAl ( J ) = SAl ( J ) / 6 . 
RETURN 



END 

SUBROUTINE MULTI (Wl , V , A, B) 

SUBROUTINE PREMULTIPLIES Wl MATRIX BY TRUNCATED A MATRIX 
(FIRST AND LAST ROWS OMITTED) , POSTMULTIPLIES PRODUCT BY 
TRUNCATED B MATRIX TfIRST AND LAST COLUMNS OMITTED) , AND 
SUBTRACTS INTERIOR ELEMENTS OF Wl FROM CORRESPONDING 
, ELEMENTS_,_OF ^Y ,, tii titii i t i ^ 

l 4 % /% /% 4 \ /C* M M X /f M A M 4 % 4 % 4 \ 4 % 4 \ 4 % 4 % 4 % 4 % 4 \ /% 4 % / 

IMPLICIT REAL* 8 (A-H ,0- Z ) 

COMMON 7 CMlA/ NLAT , NLONG 
DIMENSION Wl(ZQ) ,V(ZP) ,A(1) ,B(1) 

LATX=NLAT+ 1 
LONGM=NLONG - 1 
FCU= Wl ( 1 ) 

RCU= W 1 ( 3 “NLONG + 1 ) 

L1=4*NL0NG 
L2=L1+ 1 
L3=L1+4*LATX 
DO 2 J=2, LONGM 
K=2“ ( J- 1 ) 
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r* 



'W1(J 



(L ) + A (K+ 3 



C 

C 

520 

521 



L=3*NLONG+J 
FCC=A(K+1)*FCU+A(K 
RCC=A(K+1)*RCU+A(K 
FCU=W1 ( J ' 

rcu=wiil 

W1(J)=FC 
W1(L)=RCC 
LASTA=2*NLONG- 1 
DO 4 J= 2 , NLAT 
L=4*NLONG+J 
LL=L+LATX 
K=L+ 3 '“LATX 
KK=K-LATX 
Wl ( LL ) = A ( 3 ) *W1 ( L ) 

"■ =a(lasta}*wi(k) 

6.520) (W1(L) ,L=1, LI) 

6.521) (WliLj , L=L2 , L3 ) 

(//,' INTERMEDIATE RESULT 



A(K+3 )*W1( J+ 1) 
‘)*W1(L + 1) 



WlCKK 

WRITE 

WRITE 

FORMA 

FORMAT 



6F8 .2) 



Wl ' , / , ( 3X , 5F8 . 2 ) ) 



lastb=2*lAtx- 

NC=4*NL0NG 

Wl (NLONG+2 ) =B ( 3 )*W1 (2)+B(2) *W1 (NC+LATX+2 ) +B ( 5 ) 
l'“Wl (NC + LATX + 3 ) 

W 1 ( 2 * NLONG - 1 } = B ( 3 ) *W 1 ( NLONG - 1 ) + B ( 2 ) *W 1 ( NC + 2*LATX + 2 ) 
l+B(5)*WlCNC+2*LATX+3) 

Wl ( 2 "NLONG + 2 ) = B ( LASTB - 2 ) *W1 ( NC+ 2*LATX- 2 ) + B ( LA STB - 3 ) 
"LATX - 1 ] + B T LAS TB ) * W 1 C 3 *NLONG + 2 J 

LASTB - 2 ) ' ; W1 (NC+ 3 *LATX- 2 ) + B (LASTB - 3 ) 
+B (LASTB )*WI (NC- 1) 



l"Wl(NC+2 _ 

W1C3*NL0NG-1)=B 
1*W1(NC+3*LATX-1 
J2=NLONG-2 
DO 6 J=3,J2 

W 1 ( J + NLONG) = B ( 3 ) * W 1 ( J ) 

W 1 ( J + 2 '-NLONG) = B f LAS TB ) * W 1 ( J + 
URL=W1 (NC + LATX+ z) 

BRL= W l(NC + 2 *LATX+ 2 ) 

J2=LATX- 2 
DO 8 J= 3 , J2 
ND=2*(J-1) 

NCU=NC+LATX+J 



3 ''NLONG ) 



URC = B ( ND+ 1 ) "URL + B ( ND ) *W 1 ( NCU ) + B ( ND + 3 ) ' V W 1 ( NCU + 1 ) 
BRC = B (ND+ li“BRL+B (ND) “Wl (NCU + LATX) + B (ND + 3 ) 



C 

C 

C 



1*W1 (NCU+LATX+ 1 ) 

URL=W1 (NCU) 

BRL=W1 (NCU+LATX) 
wiTncu ) =URC 
Wl (NCU+LATX) =BRC 
W1(NC+LATX+2)=W1 (NLONG+2) 
W1(NC+2' ; LATX-1)=W1(2' V NLONG + 2 
Wl ( NC + 2 “'LATX+ 2 1 = Wl ( 2 “NLONG - 1 
Wl (NC + 3 “LATX- 1 ) =W1 ( 3‘“NLONG - 1 

CORRECT V 

NLATM=NLAT- 1 
DO 10 J= 3 , NLATM 
L=NC+LATX+ J 
K= (J-l) “NLONG+2 
V(K)=V(K)-W1(L) 

L — LATX + L 



10 



12 



C“ 

C 

C 

C 

C 



K= J -NLONG -1 
V(K>V(K)-W1(L) 
J2=NL0NG-1 
DO 12 J=2,J2 
L=NLONG+ J 

v(l)=vTl)-wi(l) 

L=2“NLONG+ J 

K= (NLAT - 1 ) “'NLONG + J 

V(k)=V(K)-W1(L) 

RETURN 

END 



SUBROUTINE FFFDB (X , E , GA , SA) 

THIS SUBROUTINE SOLVES A SUCCESSION OF ONE - DIMEN S IONSAL 
PROBLEMS. THE RELEVANT COEFFICIENT MATRIX C IS FIRST 
FORMED, THEN FACTORED, FOLLOWED BY FORWARD REDUCTION, 
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c 

c 

c 



c 

c 

c 



2 

C 

C 

C 



DIVISION BY THE DIAGONAL ELEMENTS, AND BACK SUBSTITUTION 
THE PROCESS IS CARRIED OUT NLATM TIMES. 

IMPLICIT REAL'- 8 (A-H.O-Z) 

COMMON/ CM I A / NLAT . NLONG 

DIMENSION X(l) ,E?1) ,GA(1) ,SA(1) ,C(ZU) 

NLATM = NLAT- i 
LONGM=NLONG- 1 
LONGMM=NLONG- 2 
DO 50 L=l, NLATM 

FORM COEFFICIENT* MATRIX C 

D1=E(L) 

C ( 1 ) = SA ( 2 ) •*■ D 1*GA ( 2 ) 

j2=2*nl6ng-5 
DO 2 J=2,J2 

C( J)=SA( J+2 )+Dl'"GA( J + 2 ) 



FACTOR C 




C ( 



m 



C ( 3 ) 



MP 



J2=L0NGMM-1 
DO 5 J=2,J2 
TEMP = C ( 2'- J + 1 ) / C ( 2* ( J - 1 ) ) 



7 

1000 



C(2'*'J) = C(2*J 
Itrc(2' ,; “ ' 
C(2"J+1 
GO TO 8 



. 2*J ) )7 ,7,5 
■'0+l) = TE^P 



- TEMP "'C ( 2'* J + 1 ) 



WRITE (6 ,1000)J,C( 2*J ) 
FORMAT?'/,’ STOP - MATRIX 
L’ NONPOSITIVE PIVOT FOR I 



C 

C 

C 

8 

10 

C 

C 

C 



12 

C 

C 

C 



14 

50 



2D20 . 12) 
STOP 



TRIX NOT POSITIVE DEFINITE’ ,, 
EQUATION ’,13,//,’ PIVOT = 



//,, 



PERFORM FORWARD REDUCTION 

J2= (L-l)*LONGMM 
DO 10 J=2 .LONGMM 

X(J2+J)=X(J2+J)-X(J2+J-1)*C(2*(J-1)+1) 
DIVIDE BY DIAGONAL ELEMENTS 



X( J2+1)=X(J2+1)/C(1) 

DO 12 J=2. LONGMM 
X(J2+J)=X?J2+J)/C(2*(J-1) ) 

BACK- SUBSTITUTE 

DO 14 J= 2, LONGMM 
JB= J2+LONGM- J 

XT JB ) =X? JB ) - X ( JB+ 1 ) * C ( 2“ ( LONGM- J ) + 1 ) 

CONTINUE 

RETURN 

END 
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